clear;close all;
for i = 1:1%(1:10)*1e-5
%     for d =(20:20:300)*1e-6
%          for d =(.1:.5:10)*1e-6
    
%% 输入处理
% diretory = "F:\Documents\CST\self_design_1\里程\t4cf";
%     vars = {'r12mag','r12pha','r21mag','r21pha','t12mag','t12pha','t21mag','t21pha',...
%         'epsilon1','epsilon2'};st = 0;
%     vars = {'r12re','r12im','r21re','r21im','t12re','t12im','t21re','t21im',...
%         'epsilon1','epsilon2'}; st = 1;
% diretory = "F:\Documents\CST\多重反射理论复现\多重反射理论复现 模型一 时域版\OUT DATA";
% vars = {'r12','r21','t12','t21','epsilon1','epsilon2'};st =0; % 幅值相位模式
% diretory = "F:\Documents\CST\frontier in physics\四角中心框模型的多重反射理论测试\四角中心框 频域\OUT DATA";
% diretory = "F:\Documents\CST\frontier in physics\四角中心框模型的多重反射理论测试\四角中心框 平面简化版 时域\OUT DATA";
diretory = "F:\Documents\CST\多重反射理论复现\H形复现\OUT DATA";
vars = {'r12','r21','t12','t21','epsilon1','epsilon2'};st =0; % 幅值相位模式

% st = 1;
%%
savetotal = [];
%% 循环条件也要更改
for eachturn = 1:6
    var = vars{eachturn};
    myfilename = diretory+['\',num2str(var)];%文件名
    tmp=readforinter(myfilename,eachturn);
    savetotal=[savetotal;tmp];
%{1
    %如果要显示每个变量名就打开这段程序
    makename=[var];% 先定义这样一个变量名
    eval([makename,'=','tmp',';']); 
%}
    clear tmp;
end
%     savetotal(savetotal(:,1)<0.75,:) = [];
%     savetotal(savetotal(:,1)>1.9,:) = [];
if st == 1 % re imag
    R12 = r12(:,2)+1j*(r12(:,3));
    R21 = r21(:,2)+1j*(r21(:,3));
    T12 = t12(:,2)+1j*(t12(:,3));
    T21 = t21(:,2)+1j*(t21(:,3));
else % mag phase
    M = -pi/180;
%     M = 180/pi;
%     M = 1;
    R12 = r12(:,2).*exp(-M*r12(:,3)*1j);
    R21 = r21(:,2).*exp(-M*r21(:,3)*1j);
    T12 = t12(:,2).*exp(-M*t12(:,3)*1j);
    T21 = t21(:,2).*exp(-M*t21(:,3)*1j);

end

%% 通用区
if length(epsilon1)>length(r12)
    difference = length(epsilon1)-length(r12);
    disp(['epsilon的原长度为 ',num2str(length(epsilon1))]);
    qian = floor(difference/2);
    hou = difference - qian;
    epsilon1(end-hou+1:end,:) = [];
    epsilon2(end-hou+1:end,:) = [];
    epsilon1(1:qian,:) = [];
    epsilon2(1:qian,:) = [];
    clear qian hou difference;
end

epsilon_mode = 2;
switch epsilon_mode
    case 1
        % epsilon = 3.1+3.1*.07*1j;
        epsilon = 3.42*(1+.0027*1j);
    case 2
        epsilon = epsilon1(:,2) + 1j*epsilon2(:,2);
end
% d = 250e-6;
% d = 10*1e-6;
% d = 4.2e-6;


FRE = r21(:,1);
k = 2*pi.*FRE*1e12/3e8;
BETA = sqrt(epsilon).*k*d;

R = R12- (T12.*T21.*exp(1j*2.*BETA))./...
    (1+R21.*exp(1j*2.*BETA));

figure();
subplot(411);
plot(FRE,r12(:,2));hold on ;
plot(FRE,r21(:,2));
plot(FRE,t12(:,2));
plot(FRE,t21(:,2));legend;title('mag');
subplot(412);
plot(FRE,-r12(:,3));hold on ;
plot(FRE,-r21(:,3));
plot(FRE,-t12(:,3));
plot(FRE,-t21(:,3));legend;title('phase');

% plot(FRE,abs(R))
subplot(413);

A = 1-abs(R).^2;
plot(FRE,(A));title(['吸收','d=',num2str(d)]);
subplot(414);
plot(epsilon2(:,1),epsilon2(:,2))
title('epsilon''')
end
